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'^ Pseudo Nambu-Goldstone bosons (pNGBs) are naturally light spin-zero particles, 

t^ which can be interesting dark matter (DM) candidates. We study the phenomenology 

^-^ of a pNGB 9 associated with an approximate symmetry of the neutrino seesaw sector. 

7—{ A small coupling of 9 to the Higgs boson is induced radiatively by the neutrino Yukawa 

^ couplings. By virtue of this Higgs portal interaction (i) the pNGB acquires a mass mg 

l^ proportional to the electroweak scale, and (ii) the observed DM relic density can be 

^ generated by the freeze-in of ^-particles with mass mg ~ 3 MeV. Alternatively, the 

^ coupling of 9 to heavy sterile neutrinos can account for the DM relic density, in the 
window 1 keV < mg < 3 MeV. The decays of 9 into light fermions are suppressed by 
the seesaw scale, making such pNGB sufficiently stable to play the role of DM. 
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1 Introduction 

The astrophysical and cosmological evidence for dark matter (DM) requires that a new 
particle is introduced, with constrained mass and interactions. These interactions must in 
particular maintain the DM candidate stable on cosmological time scales and generate the 
observed DM relic abundance. Many candidates can satisfy all constraints. Much fewer, 
though, provide a convincing explanation for the required values of the DM parameters. As 
examples, in most models (i) an ad-hoc symmetry is postulated to insure stability, and/or (ii) 
the DM mass is assumed to be of order (or below) the electroweak (EW) scale as demanded 
by the Weakly Interacting Massive Particle paradigm, and/or (iii) new interactions with 
strength similar to the EW ones, not motivated by other means, are just postulated. A 
scenario which is potentially predictive, and could provide motivations for some of these 
assumptions, consists in identifying the DM candidate with the pseudo-Nambu-Goldstone 
boson (pNGB) of a spontaneously broken global symmetry. 

An interesting feature of such scenario is that the DM couplings to the Standard Model 
(SM) fields, which could render the DM candidate unstable on cosmological times, are sup- 
pressed by powers of the spontaneous symmetry breaking (SSB) scale. As a result, provided 
this scale is large, the decay width of the particle can be sufficiently small. A well-known 
example is the axion, the pNGB associated with the Peccei-Quinn U{1) symmetry, broken 
spontaneously at the scale fpq. The axion is a viable DM candidate when fpq lies in the 
interval ~ (10^^ — 10^^) GeV. Another possibility is to consider SSB at the seesaw scale, 
which should be large to account for the smallness of the neutrino masses. One global sym- 
metry related to the seesaw is U{1)b~l- In the case that it is spontaneously broken by an 
EW singlet, coupled to heavy sterile (right-handed) neutrinos, the corresponding NGB is 
known as the singlet Majoron [Ij. Models with a pseudo-Majoron as DM have been stud- 
ied [2]-[I!- More generally, we will consider any (family-dependent) global symmetry whose 
spontaneous breaking contributes to sterile neutrino masses. Symmetries associated with 
the seesaw sector have the advantage that the pNGB interactions are anticipated to exist 
anyway to account for neutrino masses. 

Another remarkable property of the pNGB scenario, relevant for the DM phenomenology, 
is that the mass and the scalar potential interactions of a pNGB vanish in the limit of exact 
symmetry. Therefore, when the symmetry is explicitly broken by a unique (or dominant) 
term, the scalar interactions and the DM mass are necessarily related, and, if the relic den- 
sity is essentially determined by such interactions, this implies a one-to-one correspondence 
between the DM relic density and mass. For instance we will consider a framework where 
the source of explicit breaking generates radiatively a 'Higgs portal' interaction of the form 
X9'^WH/2, with 9 the pNGB particle and H the SM Higgs doublet. This coupling generates 
a mass ml = Xv'^ with v = 174 GeV and, at the same time, a relic density from the 9 — H 
interactions. We will show how, given the known thermal distribution of SM particles in the 
early universe thermal bath, the Higgs portal can generate a DM relic density through the 
freeze-in mechanism (see e.g. Ref. [8]). The observed relic density is obtained for a unique 
possible value of the DM mass, which turns out to be mg ~ 2.8 MeV for a Higgs boson mass 
of 120 GeV. There is also a second possibility to generate the correct relic density, by the 
freeze-in of the interaction of 9 with the sterile neutrinos. These production mechanisms 



differ from tliose previously considered for a pNGB candidate for DM. 

A tliird interesting feature of the pNGB setup is that the DM mass is related to the scale 
of explicit symmetry breaking, which can be identified with a physical mass scale already 
present in the theory, analogously to the QCD scale for the axion. For instance, in the 
scenario we will consider below, the neutrino Dirac Yukawa couplings provide the source of 
explicit breaking of the global symmetry. Therefore, the DM mass is proportional to the EW 
scale, rrig = Xv"^, and this provides a justification for the presence of a scalar DM particle at or 
below the EW scale. Moreover, the coupling A can be computed as a function of the seesaw 
couplings. In turn, we will show that in this framework the DM mass can be protected from 
quadratically divergent corrections. The mechanism we invoke to remove these corrections 
is a variation of the one proposed long ago by Hill and Ross for very light pNGBs with 
both scalar and pseudo-scalar couplings [9]. We will consider a U{l)x symmetry such that 
its explicit breaking requires several Yukawa couplings at the same time, thus lowering the 
degree of divergence of the contributions to the pNGB effective potential. This mechanism 
has similarities with the collective breaking mechanism introduced to protect the EW scale 
in little Higgs models, where the SM Higgs is the pNGB (see Ref. [lOj for a review). 

All in all, we will show that the seesaw interactions (i) are associated with global symme- 
tries broken spontaneously at the heavy neutrino mass scale, whose largeness guarantees the 
stability of the pNGB DM candidate, and (ii) have a source of explicit symmetry breaking 
built-in, and therefore induce the mass and scalar interactions of the DM. Moreover, they 
can lead to (iii) a one-to-one correspondence between the DM mass and the DM relic den- 
sity, with (iv) a justification of the presence of DM at low scale, and with (v) a DM mass 
protected from large radiative corrections. 

The paper is organized as follows. In section 2 we display the relevant effective interac- 
tions of our pNGB candidate for DM. In section 3 we show how a pNGB relic density of the 
order of the observed DM density can be obtained through freeze-in. In section 4 we present 
a class of seesaw models where the required pNGB interactions are generated naturally. In 
section 5 we compute the pNGB couplings to SM particles and analyze the corresponding 
constraints on the DM stability. We conclude in section 6. 

2 Effective interactions of the pNGB 

In this section we introduce the effective lagrangian of a pNGB 6, associated with a global 
symmetry of the neutrino sector, broken spontaneously at the seesaw scale /. We postpone 
to section |4] the detailed description of a class of models where such a light pNGB may 
emerge naturally. For the present purposes, it suffices to specify the effective interactions of 
6 at low energy (below /) . 

In the class of models under consideration, 6 couples to a sterile neutrino z/^ as follows: 

Cj, = -^z/VV^// + h.c. = -"^NN + ^9N^,N + -^O^NN + . . . , (1) 
2v^ 2 2^2 4^2/ 

where g is a. Yukawa coupling and the 4-component Majorana neutrino is defined as usual 
by A^ = {u'^ v'^'^Y . For simplicity, we consider only one sterile neutrino, whose mass ttim = 



gf /\/2 is generated by SSB. In this case the lagrangian shown in Eq. (ij) is the same as for 



the singlet Majoron model [T], where the spontaneously broken global symmetry is lepton 
number. However, we will show in section |4] that in realistic cases with two (or more) sterile 
neutrinos the relevant symmetry is not the ordinary lepton number, rather each sterile 
neutrino carries a different charge, and consequently the pNGB 9 cannot be identified with 
the Majoron. We remark that the 9 — N interaction is just an example of a NGB interaction 
with a heavy fermion charged under the associated global symmetry. As a consequence, 
the implications for dark matter phenomenology, discussed in the following sections, hold 
qualitatively in a more general set of theories, with heavy fermions other than the sterile 
neutrinos. 

A non-vanishing potential for the NGB is generated when one introduces a source of 
explicit breaking of the global symmetry. In the models we will consider, this source is 
given by a certain set of neutrino Dirac Yukawa couplings, so that the 9 effective potential 
necessarily involves the SM Higgs doublet H . The relevant, radiatively induced coupling is 
given by 

2 2 v^ 4 ' ^ ' 

where v ^ 174 GeV is the Higgs vacuum expectation value (vev), h is the physical Higgs 
boson and we adopted the unitary gauge. Note that the mass of the pNGB 9 is proportional 
to the EW vev, rrig = Xv"^. This relation is sometimes referred to as the 'conformal' limit, 
since 9 has no bare mass term, rather the Higgs portal interaction A generates both mo and 
the 9 — h couplings. We remark that the phenomenology (e.g. the DM mass, its relic density, 
etc.) induced by the Higgs portal in Eq. ^ does not depend on the details of the underlying 
physics which generates radiatively A. 

In the specific seesaw models that we will build, A is generated by logarithmically diver- 
gent neutrino loops, involving both the 9 — N coupling g as well as neutrino Dirac Yukawa 
couplings, denoted generically with y. It turns out that the relation between the mass of 9 
and the neutrino parameters can be written schematically as 

2 X 2 2 2 2 log(AVm^) 
mQ = \v ^g y V —- , 3 

where A is some cutoff scale at or above /. The light neutrino mass scale is given by the 
standard seesaw relation, m,^ = y^v^ jm^. It is useful to use Eq. ([3| to express the coupling 
(7 as a function of the relevant energy scales in the theory: 
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We will take the factor k = log(A^/m^)/(87r^) to be of order one. Then, for a fixed value 
of m.^, r/iTv and also / = \/2m,N/g are determined as a function of m.g and g, which are 
basically the only two free parameters. We require the scale / to be smaller than the Planck 
scale, Mp = 1.22 x 10^^ GeV. In addition, we assume for simplicity that / is larger than 
one TeV, in order for the effective theory to contain only the SM particles, the pNGB 
and the sterile neutrino, while the degrees of freedom involved in the SSB are decoupled. 



These two requirements exclude the green shaded regions in the {rrig — g) plane shown in 
Fig. [ij The dotted lines in Fig. [l] correspond to constant values of the sterile neutrino mass, 
ruN = ml/{g^m^k) = 10^ 10^ 10^°, 10^^ GeV, from top to bottom. 

3 Relic density from freeze-in of the pNGB 

We now study the role of the interactions in Eqs. ([I]) and ^ for the production of a 6* relic 
density, which can play the role of DM. In this section, we implicitly assume that 6 is stable 
on cosmological time scales. The region of parameters where this stability is achieved will 
be analyzed in detail in section |5j There are several mechanisms to produce the DM relic 
density. We begin by briefly recalling these various possibilities and confront them with the 
parameter space of our scenario. 

A well-known way to produce a cosmological density of a DM particle species is the usual 
freeze-out: it requires an annihilation cross section large enough to overcome the expansion 
of the universe and thermalize the DM particle. The interactions of 6 come from the terms 
in Eq. ([I]), which we refer to as the 'sterile neutrino portal', and the terms in Eq. (pi), that 
is, the 'Higgs portal'. These interactions may or may not lead to the thermalization of 6 in 
the early universe, depending on whether the interaction rate F gets larger than the Hubble 
rate H. If it does, later on 6 necessarily decouples from the plasma and, if stable, acts as 
DM with a certain energy density p^, which depends in general on its mass and annihilation 
cross section. 

Before determining pg numerically, it is instructive to identify the thermalization region 
in first approximation, by evaluating the 6 interaction rate F and requiring T > H. We 
start with the Higgs portal, considering for simplicity the h ^ 99 decay only, because in our 
scenario decays turn out to dominate with respect to scattering processes (see Fig. [2] below). 
The decay rate is 

nh^99) = ^X^^Jl-^^. (5) 

167r m/j y m^ 

This should be compared to the Hubble parameter in the radiation epoch, which is given 
hy H = 1.66a/^ T^/Mp, with gP = 106.75 the degrees of freedom corresponding to the SM 
content (we shall assume such value in the rest of the paper). Taking T ~ irih to be the 
relevant temperature, we find that, for a coupling 

the decays and inverse-decays are able to produce a thermal population of 6'-particles. In 
other words, given the m?Q = Xv"^ relation, one needs nig > 44 MeV, which translates into 
the thin vertical line in Fig. [ip] 



^Notice that Eq. ([6| is the condition for DM to thermahze with the SM sector at T ~ rrih (where the rates 
reach the maximum, see Fig. ^ below) . This condition is different from the one usually considered for the 
freeze-out scenario, where one requires that the scattering processes are in thermal equilibrium at T < mg. 
The latter condition leads to a value of A about two orders of magnitude larger than the one in Eq. (l6|) . 



Similarly for the interaction with the sterile neutrino A^, the process NN — )■ 66 has a 
rate 

4 

r = {aVrei)nN ^ ;73T ^ ^2x 10"^ g^rriN , (7) 

256tc mjj 

with riAT the sterile neutrino number density, and we used T ~ ttitv- This leads to thermal 
production of 6'-particles provided 

eV \ ^/^ f me \ Vs 
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where we made use of Eq. Q. This gives the thin diagonal hne in Fig. [I] We do not discuss 
the effects of the couplings of 6 to light fermions, because in the models under consideration 
they turn out to be far too small to thermalize 6 (see section [s]). 

Outside the thermalization region, where F is smaller than H, one may still have a final 
Qg matching Qdm (as usual we define Qi = Pi/pcrit, i-e. the ratio between the energy density 
Pi and the critical density of the universe). The idea is that in the early universe a population 
of ^-particles is produced through pair production scattering processes, but with a density 
smaller than the thermal one. For appropriate values of tuq one may still reach a relic density 
equal to floM- This process, called freeze- in, actually happens in our scenario as we shall 
describe below in detail (for generic properties of the freeze- in mechanism see Ref. [8]). 

In the rest of this section we first introduce the ingredients to make a numerical calculation 
of the 6 relic density, and we then apply them to the freeze-out and freeze-in processes. 
Afterwards we also discuss some non-thermal production mechanisms and end up with a 
summary of our results. 

3.1 Boltzmann equation and reaction densities 

From the Higgs portal, ^-particles are created from the following scattering processesjj 

a{hh^66) : 

a{WW^66) -- 

a{ZZ^66) -- 

a{ff^66) - 

where Th is the total Higgs decay width and tIc is the number of colours of the fermion /. 



A2 f s - Am^y^' f s + 2miy 


(9) 


327rs V'^ ~ 4??i^/ V -^ ~ ""^ft, / 


1 A^ fs-AmlY^ s^ - Asml^ + 12m^^ 


(10) 


932ns \s- Ami,) (^ - K)' + ^^H ' 


1 A2 /s-AmjY^ s2-4sm| + 12m| 


(11) 


9 327rs Vs - 4m|y (s - ml^ + mlVl ' 


1 1 A^ m){s-Am)Y/\s-Amlf/^ 


(12) 


4nc levrs {s - mlf + mlVl 



^In Eqs. (|9])-( 13 ) we included the usual initial-state spin and color average factor, as well as the final state 
identical particle factor 1/2. In a{hh — > 69) we neglected subleading terms of order A'^ and A** (which are 
displayed e.g. in Ref. [H]). 
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Figure 1: The pNGB coupling to the sterile neutrino, g, versus the pNGB mass, mg, in GeV. We 
fixed rriu = 0.05 eV and k = 1. The upper (lower) green shaded region is excluded by requiring 
the SSB scale / to be larger than one TeV (smaller than the Planck scale). The four dotted lines 
correspond to constant values of the sterile neutrino mass: from top to bottom, m^ = 10^, 10®, 
10^*^, 10^^ GeV. The thin diagonal (vertical) line indicates the lower value of g (mg) that leads 
to a thermalization of 9. The thick line corresponds to a ^ relic density equal to the observed 
DM relic density, as follows from the numerical solution of the Boltzmann equation (for ruh = 120 
GeV). Since 6 is produced by freeze-in, its relic density grows with its couplings g and A = uIq/v'^, 
therefore the region below and to the left (above and to the right) of the thick line corresponds to 

^9 < ^Dm{> ^Dm)- 



Similarly, ^-particles can be created through the scattering with a sterile neutrino, 
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where we neglected mg in front of m^v and ^/s and we defined /Sat = (1 

In the following it will be convenient to consider separately the effect of the Higgs de- 
cay, Eq. ([5]), which means that, in order to avoid double counting, the on-shell part of 
the scattering processes must be subtracted (whenever there is one). In the narrow width 



approximation this means to perform in Eqs. (10)-(12) the substitution 



TUhT 



h 



mlY + m?r? 



-)■ 



rrihTh 



Tr6{s — ml)6{s — 4m^ 



(14) 



W, Z, f. In the following, all scattering quantities we will consider are meant to be 



for i 

the subtracted ones, using Eq. (|14|). To calculate the relic density that one obtains from 



these processes, either through freeze- in or freeze-out, one has to integrate the Bohzmann 
equation. 



zH(z)s(z)Y'(z) 



1 



Yoiz' ^ ' 



rr(^l 



ilDiAz) + lannih{z)) , (15) 



with z = rrih/T conventionally taken as the evolution parameter. Here H{z) is the Hubble 
parameter, s{z) the entropy density, and Yq = ng/s with uq the number density. The 
reaction density 7^^ contains the effect of the Higgs boson decay, while ■jannih is the sum of 
the (subtracted) reaction densities of the scattering processes in Eqs. ([9])-(13). The reaction 
densities are given by 



Nenl'^.nh^ee), (16) 



ID^ 



7(a 6 o 1 2) = Ne J J dpjp.rjf^" J J dp,dp2{27r)^6\pa + Pb - Pi - P2) \M\^ 

involving the Bessel functions Ki^2- We have defined dp = d^p/{{2TT)^2E). Here A^g = 2 is 
the number of 9 particles produced per decay or annihilation, f^'' = (e^'/^ ± 1)^^ ^ e^^'/"^ 
is the Maxwell-Boltzmann energy distribution, |A^p is the amplitude squared summed over 
initial and final spins (with no averaging), Smin = max[(?7ia + rrife)^, {mi + ^12)"^], and the 
reduced cross section is defined by 

a(a6 O 1 2) = ^ ^[is-<-<?-^rnl<] ,(,,^12), (18) 

Cab S 

with a the particle physics cross section of Eqs. (|9])-(13), ga,b the number of degrees of 
freedom of the particles a,b and Cab a combinatorial factor equal to 2 (1) if a and b are 
identical (different). 

3.2 Freeze-out 

As it is well-known, if 6 freezes out relativistically, that is at T > mg, its relic density is 
independent of the annihilation cross section, because it decouples when the thermal number 
density, Y^'^ = n^^ / s is still independent of the temperature. The formula for the relative 
density VLq is simply 

n„,4=^78i|^. (19) 

with gl the number of degrees of freedom contributing to entropy at 9 decoupling. This 
should match the value 

Vtouhl = Q. II ±Qm . (20) 



In these equations, ho is the reduced Hubble constant Hq = Ho/{100 Mpc km s^^) ^ 0.70. 
As a result, the observed relic density can be obtained only for one DM mass value: mg ~ 
0.15 keV, where we took g^ = 106.75, valid for a decoupling temperature above 100 GeV. 
Such a light 6 can be only thermalized by the interaction with sterile neutrinos, which 
decouple at a temperature just below mjy, leading to relativistic 6''s. This result corresponds 
to the left vertical branch of the thick curve in Fig. [T| which is obtained by integrating 
numerically the Boltzmann equation. 

If instead 6 is heavier and thermalizes, to have a small enough relic density requires that 
it decouples non-relativistically. This can occur through the Higgs portal interaction. In this 
case the relic density essentially only depends on the annihilation cross section which must 
take the value (aVrei) — 10~^^ cm^s~^. To get a sufficiently small relic density, one needs to 
go to much higher values of mg, where the coupling A = mg/v"^ becomes large enough. The 
observed relic density is obtained for a unique value of nig, because of the relation between 
mg and A, which is a consequence of the pNGB nature of 6. The freeze-out through the 
Higgs portal was considered in Ref. [12], assuming this restrictive relation, and it was found 
that the DM relic density is obtained with m^ = 50 - 70 GeV (for m/, = 120 - 180 GeV). 

However, as we will see in section [5} the freeze-out values, mig = 0.15 keV as well as 
m,g = 50 — 70 GeV, require values of parameters which are excluded by the instability of 6 
on cosmological time scales. Therefore in our scenario both relativistic and non-relativistic 
freeze-out is not a viable way to produce DM. Instead, the freeze-in process turns out to be 
efficient in an allowed region of the parameter space, as we now discuss. 

3.3 Freeze-in 

In the freeze-in process 6 is produced by the annihilation or decay of a heavier particle X, 
with rates small enough not to thermalize, leading to a less-than-thermal DM population, 
which reaches a plateau at T ~ mix because at smaller T the number of X-particles is 
Boltzmann suppressed. If the X-particles are the SM particles, as it is the case through the 
Higgs portal, their number densities are known as a function of the temperature (simply 
given by their thermal distribution down to a temperature well below their mass). As a 
consequence, the created relic density depends only on the magnitude of the portaljj 

Let us ffist consider a value of the coupling g sufficiently small to neglect 6 production 
at the seesaw scale (we will see that this is the case for g < 10^'^ — 10""^, depending on m^). 
In this case only the Higgs portal interactions can account for the observed relic density. 
As said above, this mechanism is extremely predictive since, given the 'conformal' relation 
mg = Af^, the rates and hence the relic density depend only on the parameter mg. 



Integrating numerically Eq. (15), it turns out that the observed relic density in Eq. (20) 
is obtained for 

mg = (2.76 - 2.86) MeV , (21) 

where we took ruh = 120 GeV for definiteness. This result corresponds to the right vertical 
branch of the thick curve in Fig. [1} 



^For production of DM through the Higgs portal along the freeze-in scenario, see also Ref. 
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Figure 2: The decay therinalization rate, "fDh/{n^^H) (black), compared to the various scattering 
thermalization rates, lannih/i''^V ^) ^°^ '^ ~ W, Z,h,t,b (in red, blue, green, orange and purple 
respectively), as a function of z = rrih/T and for nih = 120 GeV, mg = 2.8 MeV. For freeze-in these 
rates remain always well below one, which is the thermalization value. 



An important difference between freeze- in and freeze-out is that, while the latter is gen- 
erally dominated by the off-shell annihilation, the former is naturally dominated by the 
on-shell annihilation, i.e. by the decay of the mediator involved in the annihilation. For a 
freeze-out, unless the mediator has a mass ~ 2mDM, the decay has a small effect because, 
when the off-shell annihilation goes out-of equilibrium, it has already decoupled (being more 
Boltzmann suppressed). For a freeze- in instead, the DM production occurs mostly at T of 
order the mediator mass, or slightly below, where the Boltzmann suppression is still mild. As 
a result the decay, which involves less (small) couplings, is naturally dominant, 7d^ > 'jannih- 
This is shown in Fig. [2] where the various reaction densities are plotted for rrih = 120 GeV 
and mg = 2.8 MeV. In the approximation where only the decay is included, the Boltzmann 
equation can be integrated analytically, giving [S] 



y^^AT, 
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(22) 



where qh 
1%. 



1. This analytic result turns out to differ from the numerical result by less than 



To understand Eq. (22) let us first discuss why the freeze-in is infrared dominated, with 

< 



most of the 6 production occurring at T < nih (just before rf^^ gets too much Boltzmann 
suppressed). The 6 number density produced at temperatute T is essentially given by the 
number of decays per unit time occurring at this temperature, rfj!^{mh/T)T{h — )■ 99), times 
the number of ^-particles produced per decay, Ng, times the Hubble time, 1/H . Therefore, 
taking into account thermal averaging properly, one has Yg ~ 'y£)h/{Hs), which is nothing 



10 



but the rate of thermalization 'jj:)h/{n'g'H) shown in Fig. pi times the relativistic value of 
F/^. Thus one obtains Yg ex MpmhT{h — )• 99)/T^, which increases as T decreases. 

Hence one finds Yg = c[Ngn'i^r{h — )• 66)/{Hs)]T=mhJ with c a numerical factor which turns 
out to be equal to {37r/2)/K2{l) — 2.9. It is larger than one because, as Fig. |2] shows, the 
maximum value of the decay reaction density occurs at T ~ m/i/3.5 rather than at T ~ mhrj 
With this in mind, and taking into account the dependence T{h — )■ 99) oc X'^/rrih oc mg/rrih, 
the DM relic density from the Higgs decay scales as ml/m^, so that the value of me needed 
varies very little with the exact value of the relic density, and the larger the Higgs mass is, 
the larger mg has to be. For example for rrih = 140, 180, 300 GeV one gets the observed relic 



density for mg = 3.0, 3.6, 4.9 MeV instead of 2.8 MeV in Eq. i^. 

It is interesting to discuss what is the contribution of each scattering channel separately, 
even though they have a small effect with respect to the decay. As can be seen from Fig. [2| 
among the scattering channels, the W one gives a contribution slightly larger than the Z 
and h ones, and much larger than the b and t ones. Numerically, for rrih = 120 GeV, 
the W, Z, h, b, t scatterings contribute to the total number of ^-particles produced in the 
proportions 1 : 0.4 : 0.4 : 4 x 10~^ : 10~^ respectively, while the W scattering alone gives a 
contribution 500 times smaller than the decay channel. 

Similarly to the decay, at a temperature T > nta., a scattering process aa — t- 99 gives a 
number density Yq ~ lannih/i^^)y which is infrared dominated too. Again, this is nothing 
but the thermalization rate lannih/i^^H)-! displayed in Fig. |2| up to the multiplicative con- 
stant Yq'^ . One gets Yq = Ca[yannih/ (^ ^)]T=ma with Ca & numerical factor which, unlike for 
the decay, does not take a unique value. It depends on nig as well as on the s dependence 
of the scattering cross section (which fixes the position of the peak of the reaction densities 
in Fig. [2]). For the W, Z and h channels Ca — 2.2, whereas for the fermion channels q ~ 1.7 
and Cb — 6. One understands consequently why the W, Z and h channels, which give con- 
tributions Yg oc Mp/my/^z,hi dominate over the b channel, which gives Yg oc mlMp/m\. As 
for the top channel, it is suppressed too because for s > Ami i^ i^ more threshold suppressed 
than the other channels, see Eqs. ([9])-(12), and for higher s it gets quickly suppressed by the 
ml/ s^ asymptotic behaviour of the top cross section (as compared to the 1/s asymptotic 
behaviour of the W , Z and h cross sections). 

Let us now move to the region mg < 3 MeV, where the Higgs portal alone would lead 
to a too small relic density. In this case the N annihilation process can do the job, as long 
as the reheating temperature is large enough to produce a thermal population of sterile 
neutrinos, which we assumejj The annihilation rate depends on two free parameters only, 
g and ttiq, with the sterile neutrino mass ttijv determined by Eq. (|4]). From the freeze- in of 



^Actually if one takes the value of ^Dh/{Hs) at its maximum value, i.e. at T ^ mii/3.5, one obtains a 
number density equal to 70% of the exact result. Therefore, the production of a particle through freeze-in of 
a decay is essentially given by the maximum value of the corresponding thermalization rate. Note also that 
if one takes instead Ngnff^T(h — > 99)/{Hs) at its maximum value, which occurs for T ^ mh/3.0, it turns out 
that one gets an even better approximation, as this overestimates the exact result by only 3%. 

^It is well-known that, if a sterile neutrino has Yukawa couplings inducing a light neutrino mass of the 
order of the atmospheric or solar neutrino mass scales, the decays/inverse decays of the sterile neutrino are 
in thermal equilibrium at T ~ rriN- There are various cosmological issues which can put constraints on 
thermal sterile neutrinos, see e.g. [2]. However they become relevant for sterile neutrino masses below the 
~ GeV scale, which are not pertinent in our case, see Figs. 3J4 below. 
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the reaction NN — )■ 99 one obtains as usual Yg ~ ['yannih/{Hs)]T=mN- ^^^ then finds that 
the DM rehc density scales as g^m^/me, as confirmed by the numerical integration of the 
Boltzmann equation. 

Our numerical result in the {mg — g) plane is shown in Fig. [I] for rrii, = 0.05 eV. The 
diagonal branch of the black thick line corresponds to a freeze-in through the sterile neutrino 
portal, with the correct value of the DM relic density. The required value of the 9—N coupling 
is given approximatively by (? ~ 2- 10~^ (mg /MeVY^^ {eV / Tn,^y^^ . Above (below) the line the 
relic density is too large (small). Since the relic density is proportional to the sixth power 
of g, the predicted value of g depends very weakly on the cutoff A, which we fixed in all the 
numerical analysis, k = log(A^/'m^)/(87r^) = 1. 

3.4 Non-thermal production 

In general, one may think of non-thermal, more model dependent, 9 production mechanisms. 
A class of mechanisms is related to the phase transition at the scale of SSB /. At this epoch 
(part of) the energy stored in the scalar potential false vacuum may be transferred to the 
pNGBs. There could also be cosmic strings produced at the phase transition, which would 
decay into pNGBs, but their exact density depends on model details. Here we assume that 
these contributions to Qg are negligible. This is the case, in particular, if the universe 
underwent an infiationary phase at temperatures smaller than / (note that the reheating 
temperature can still be larger than uin = gf 1^2^ since the coupling g is much smaller than 
one in our scenario). 

A potential source of non-thermal production of pNGB dark matter could be provided by 
the oscillations of the field 9 around the minimum of its effective potential, in case the value 
of 9 at high temperature is displaced from such minimum, a mechanism well-studied in the 
case of the axion. Let us summarize the axion case. The axion mass is very suppressed above 
Kqcdi therefore at high temperatures the axion behaves as an exact NOB. In this case, at 
the end of the associated Peccei-Quinn phase transition at scale /pg, the value of the axion 
field a can lie in any of the equivalent vacua described by < a/ fpq < 27r. Later, when 
ma{T) becomes important with respect to the Hubble parameter H[T), the field a begins 
to oscillate around zero, producing a coherent state of particles at rest with an associated 
Qosaii ^£qj, ^ review see Ref. [11]). One may naively think that this picture applies also to our 
scenario, with a negligible mg{T) at high temperature (where the EW symmetry is restored), 
and the oscillations beginning only at T ~ TeV, when nig is generated. If this were the case, 
taking the values of / and mg in the range we considered for the freeze-in production, one 
would conclude that f2^*^*" overdoses the universe, unless the initial value of 9 is tuned to 
be much smaller than /, or infiation takes place at temperatures below the TeV scale. 

However, the temperature dependence of m,g{T) is very different from the one of the 
axion, whose mass is generated non-perturbatevely by an anomaly at Aqcd- In our case, mg 
is generated, instead, by an explicit breaking of U{l)x and therefore it does not vanish at 
high temperatures. On the contrary, we expect m,g{T) to receive large thermal corrections. 
Even if the EW symmetry is restored for T > TeV, there are contributions to m,g that are 
not proportional to the Higgs vev v, such as 6m,g ~ (7^y^A^/(167r^)^ (see the discussion at the 
end of section E|. These lead to m^^iT) ~ constant ■ T^ at high temperature. A quantitative 
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estimate of a possible non-thermal production of 9 would then require a non-trivial study 
of the thermal evolution of the field 9 during and after the phase transition at scale /. 
Still, we just notice that mg{T) is typically larger than H{T) already at T ~ /. Thus, we 
argue that the field 9 does not acquire a random initial value of order /, but rather it sits 
in the minimum of the potential already at high temperature, or in other words coherent 
oscillations are strongly suppressed. On the basis of this argument, we assume that fi^'^'^*" is 
negligible, and that the thermal freeze-in dominates the DM relic density. 

3.5 Summary 

Our results on the 9 relic density are summarized in Fig. [T| where we show the line Qg = ^dm 
in the {mg — g) plane. The behaviour of the line with the correct relic density can be easily 
understood. For nig ~ 0.15 keV, g must be larger than (or equal to) the value needed 
to thermalize 9. For progressively larger mg, the value of g should correspond to less and 
less thermalization. For 0.15 keV < mg < 3 MeV, one can see that the required coupling 
g is progressively smaller than the thermalization value, indicated by the thin diagonal 
line. When one approaches nig ~ 3 MeV, g becomes rapidly smaller as the Higgs portal 
begins to contribute to the 9 number density. Once the Higgs portal produces by itself the 
correct amount of 6'-particles, g must be small enough to make the sterile neutrino freeze-in 
negligible. 

We conclude that the observed DM relic density can be obtained thermally for 

0.15 keV < me < 3 MeV , (23) 



where the exact upper value depends on the Higgs mass, as discussed in section |3.3[ Note 
that this upper bound holds as an absolute prediction as soon as the Higgs portal dominates 
the DM production. This is necessarily the case, in particular, if the reheating temperature 
lies below the sterile neutrino mass scale. More generally, the prediction mg ~ 3 MeV holds 
for any pNGB whose Higgs portal interaction A gives the dominant contribution to its mass 
and relic density, independently of the associated global symmetry and of the source of 
explicit symmetry breaking which induces A. In particular the SSB scale may be different 
from the seesaw scale. 

Of course the prediction for mg would change if the number of DM particles was not 
negligible already at temperatures higher than the EW scale. In this case one has to produce 
less of them through freeze-in, which means that mg has to be smaller. In other words, 
mg ~ 3 MeV constitutes an absolute upper bound for the freeze-in mechanism, and it holds 
as soon as the freeze-in production dominates over the initial population. 

4 Approximate symmetries of the seesaw sector 

In this section we consider the SM augmented with sterile neutrinos u'ji, with a global sym- 
metry broken spontaneously at the seesaw scale, and we discuss in some detail the generation 
of a mass for the associated pNGB. 
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Let us consider the most general Yukawa interactions to be added to the SM in the 
presence of gauge singlet fermions, 

- £,c = l^ruajiy^ ( 7 ) + 2"'^'^"^ + h.c. , (24) 

where la are the lepton doublets and H is the Higgs doublet, whose neutral component 
acquires a vev w = 174 GeV. Here and in the following, the mass parameters rriaj and Mij 
are intended in a generalized sense as dynamical scalar fields that may or may not acquire a 
non-zero vev. Thus, the above lagrangian has a global U{1)l symmetry with the following 
lepton number assignments: 

L{1) = +1 , L(z/") = -1 , L{H) = L{m) = , L{M) = +2 , (25) 

where we dropped flavour indices. When M acquires a vev, U{1)l is spontaneously broken 
and a massless NGB appears in the spectrum of the theory, the singlet Majoron ^1]. One 
can write more explicitly 

M,, = (7.,$ , $ = -^e^'^f , (26) 

where $ is a complex scalar with L = 2, the vev (p) = / breaks spontaneously lepton 
number, and 6 is the Majoron. 

There are various possible sources of [/(l)^ explicit breaking coming from other sectors of 
the lagrangian. There may be soft terms in the scalar potential involving $ which break lep- 
ton number; their mass scale, whose size is arbitrary, determines the induced Majoron mass 
(see e.g. f7|). Another possible source of U{1)l explicit breaking are quantum gravity effects 
at the Planck scale Mp, which can break in general all non-gauge symmetries. Assuming 
these effects are suppressed by powers of Mp, they can be used to generate a Majoron mass 
at the keV scale [21 13], which has been extensively studied as DM candidate [31 HI |5|, [6] . 

In this paper we will consider a different source of explicit breaking of global symmetries, 
provided by the set of the Yukawa couplings [9J . This allows to relate the size of the pNGB 
parameters to the fermion mass scales already present in the theory. We will focus on 



symmetries of the lagrangian in Eq. (24) other than U{1)l, with each lepton carrying in 
general a different charge. In this case the symmetries are respected only by some matrix 
elements maj and Mij, and they are explicitly broken if some other matrix elements are non- 
zero. The mass and couplings of the pNGBs will be completely determined by the seesaw 
parameters and by the choice of the cutoff, since they can receive cut-off dependent quantum 
corrections. 

The pNGB mass is, in general, quadratically sensitive to the cut-off. In order to under- 
stand the origin of quadratic divergences and the mechanism to remove them, consider first 
the explicit breaking of a U{1) symmetry in a theory with only one sterile neutrino u'^: 

-£,c = lz/^(M,e^^/^+M,)i/^+h.c. = ^{Ma+M,)NN-'-^eNj,N-^9mN+0{i9/ff) , 

(27) 
where Ma,b are real mass parameters. Here Ma is generated by the spontaneous breaking of 
the U{1) symmetry associated with the NGB 6, while Mb breaks this symmetry explicitly. It 
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is instructive to compute the fermion loops generating the pNGB mass term, m^ 6/2, using 



the right-hand side of Eq. (27): if Mb is zero, the two relevant one-loop diagrams cancel each 
other, as expected for an exact NGB. However, in the presence of the explicit breaking, a 
non-zero quadratically divergent contribution is left, which is given by 

"^e~^-^A . (28) 

The effective theory below the scale / may still contain a light scalar 6', since nig can be 
parametrically small, but its NGB nature is obscured by the quadratic dependence on the 
details of the ultraviolet completion. 

In the presence of more than one fermion family one can define several family-dependent 
f/(l) symmetries. In general, they are explicitly broken by some mass matrix entries. How- 
ever, certain ?7(l)'s will be broken only when several entries are non-zero at the same time: 
this is the key to reduce the degree of divergence of the radiative contribution to m^, where 
^'s are the associated pNGBs [9]. This fact can be understood in the language of the ef- 
fective potential K//, considering the full (active + sterile) neutrino mass matrix M.: the 
term in V^ff quadratic in the mass matrix, ~ Tr(A^A^"'')A^, is invariant under certain U{1) 
symmetries, that is, it does not depend on the associated pNGBs. Thus, the potential of 
these pNGBs contains at most terms quartic in the mass matrix, ~ Ti{A4.M.'^ M.M.'^) log A^, 
which are only logarithmically sensitive to the cutoff (for an application of this idea to eV 
scale sterile neutrinos, with the pNGB playing the role of dark energy, see [16]). Of course 
non-abelian symmetries are also possible, with several NGBs and, potentially, qualitatively 
different phenomena, but this extension will not be needed for our purposes and will not be 
considered in this paper. 

In order to identify the combination of matrix entries that induces a pNGB mass, we 



rewrite Eq. (24) with the replacement la{H/v) — ?■ Vai in the minimal case of two sterile 



neutrinos, which are sufficient for realistic light neutrino masses: 

- ^- - "' ('»- ™-) ( :n + >' "» ( t iz ) ( ;0 + ''■^- • p^' 

Here a sum over n/ = 3 active flavours is understood [a = e, fi, r). Let us identify the matrix 
entries whose phases can be removed. There are 2nf + 3 mass terms, which in general have a 
phase. One can absorb n/ -t- 2 of these phases by redefining the 'nf + 2 neutrino fields present 



in Eq. (29). As a result, there remain nj + 1 complex matrix entries in Eq. (29) (either nj 
phases in m and one in M, or equivalently nj — 1 phases in m and two in M). When one 
sets these Uf + 1 entries to zero, there are no physical phases left. 

Now, suppose these Uf + 1 complex entries vanish and the non-zero entries are generated 
by a set of scalar fields $a = p^ x exp{i9a/ fa)/V2, acquiring a vev {pa) = fa- Then, by an 
appropriate redefinition of the lepton fields, one can remove all the phase fields 6a from the 



Yukawa lagrangian of Eq. (29). Thus, these fields will have only derivative couplings to the 
leptons, resulting from the redefined lepton kinetic terms. This means that the lagrangian 
has exact U{1) symmetries, broken spontaneously by ($)'s, and 6a are exact NGBs. 

When one of the zero matrix entries is switched on, however, it is no longer possible to 
remove all the phases. This means that one ^(1) is explicitly broken, with the associated 
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pNGB 9 acquiring non-derivative couplings. Thus, a non-zero irtg is generated by neutrino 
one- loop diagrams. To estimate mg, consider the set of the (nj + 3) non-zero matrix entries. 



One can check that 9 could be rotated away from Eq. (29) if any out of four of these entries 
were switched off. This means that a f/(l) symmetry is recovered when any of them is put 
to zero, and therefore itiq must be proportional to the product of the four entries. As a 
consequence, quadratically divergent contributions to uiq turn out to be absent. 

The 9 mass can be controlled by the entries of M only, by those of m only, or by both. 
Let us specify three models with a U{l)x symmetry, which are representative of these three 
generic possibilities: 

(i) The symmetry is broken explicitly in the singlet neutrino sector. Two independent 
entries of the Majorana mass matrix M are allowed by f/(l)x, while the third is 
forbidden. For example take X{yl) = 0, ^(z/g) = 1 and one scalar field $ with 
X{^) = —2. Then Mn is allowed and M22 is generated by the vev of $, while 
M12 = 0. The NGB coupling to neutrinos reads M22 Q'^p{i9 / f)i'2^2/'^- (The Dirac 
neutrino sector is not relevant here: the m entries may or may not respect U{l)x, in 
any case the leading symmetry breaking effects are controlled by M.) When M12 is 
different from zero, neutrino loops give a non-zero contribution to mg as long as they 
contain four appropriate mass insertions: 



1 MUM12M22M12 A^ 



"^e ~ ^T 79 log^> (30) 



where we included a loop suppression factor and fi is the renormalization scale, which 
can be taken of the order of the sterile neutrino masses. The cancellation of quadratic 
divergences as well as the structure of the non- vanishing contribution to mg are dictated 
by our symmetry argument, and they can be both explicitly checked by computing the 
loops. 

The symmetry is broken explicitly within the Dirac neutrino sector only. The entries 
on a given row of the Dirac mass matrix m are both allowed, while in the other rows 
only one entry is allowed. For example take X(z/e) = 1, X(z/^^^) = 3, X(z/f) = — 1, 
X(z/2) = 1 and a scalar field with X(0) = —2. Then 171^1 is allowed, 771^2, "m^i and 
niT-i are generated by the vev of 0, and m^2 = '^t2 = 0. (For simplicity we suppose 
that the scalar fields and $ that can contribute to m and to M belong to different 
sets. The singlet neutrino sector is assumed to respect U{l)x'- in the present example 
M12 is allowed, thus giving an equal mass to the two sterile neutrinos.) This case 
reproduces closely the scenario discussed by Hill and Ross [0], which dealt with the 
quark Dirac mass matrix. Once the entries explicitly breaking U{l)x are switched on 
in m, one obtains 



") 



2 1 maim^2m^2mpi A^ 



m^-^ "^ ,7^^ "M og^, a^/3. (31) 

(iii) The symmetry is broken explicitly only by the interplay of the Dirac and Majorana 
neutrino mass matrices. One possibility is that U{l)x allows for all the entries of 
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M, e.g. taking the charges X(z/^) = —1, ^(i^g) = 1 ^^d X($) = 2; this symmetry 
can allow only one between niai and ma2, depending on the charge of //„. Another 
possibility is that f/(l)x allows for all the entries of m, but it forbids two entries of 
M, e.g. taking the charges X(z/a) = X^uf) = and X(z/|) = — X(0) 7^ 0. As usual 
the mass of the NGB 6 is generated when one zero entry is switched on. One finds 



2 1 MijMkimaima2 , A^ 
" ~ log — 

87r2 /2 ^ /i2 



^e ~ ^^^ 7^ log ^ . tj j^kl . (32) 



In order to roughly quantify the energy scales under discussion, let us indicate with M = gf 
{m = yv) a generic seesaw (electroweak) scale, possibly suppressed by a small Yukawa 
coupling g (y) with respect to the U{l)x SSB scale / (the EW scale v), and let us introduce 
also the neutrino mass scale m,y = m'^/M. Up to a loop suppression, in the three cases mg 
is of order (i) {M/ f)M = gM, at or below the seesaw scale; (ii) {m/ f)m = grrty, at or 
below the neutrino mass scale; (iii) {M/ f)m = gm, at or below the electroweak scale. In 
the following we will concentrate on the last possibility, since we are interested in scalar DM 
candidates below the electroweak scale. 

Before moving to the other properties of the pNGB relevant for DM, some comments 
are in order to assess the soundness of the above estimate for the pNGB mass me. We 
have seen that, in general, in the presence of explicit breaking of the associated global 
symmetry, mg can be sensitive quadratically to the cutoff of the theory. However, if the 
global symmetry is broken only by the contemporary presence of several couplings, then 
m,g shall be proportional to the product of all these couplings. In this case the Feynman 
diagrams contributing to the pNGB mass will have a lower degree of divergence and thus 
quadratic divergences vanish. This possibility was often employed in model-building in the 
past, and recently was extensively used in the context of little Higgs models, under the 
name of "collective breaking" (see Ref. [10] for reviews). The motivation is to stabilize the 
electroweak scale against large quantum corrections. In these models the Higgs is a NGB of a 
symmetry broken spontaneously at the scale few ~ TeV, whose mass is generated at one-loop 
level by two couplings ?/i2, whose contemporary presence breaks explicitly the symmetry. 
As a consequence, m,\ ~ |/^?/|/g^log(A^/yU^)/(167r^), that is, the sensitivity to the cut-off A 
is only logarithmic. However, in general, the quadratic divergence reappears at the two-loop 
level, with a correction to the Higgs mass Sm^l ~ y^y|A^/(167r^)^, which indicates that A 
cannot be larger than ~ Aiifew, in order for the theory to remain natural. 

It is instructive to compare this little Higgs scenario with our scenario, where the DM 
candidate 6^ is a pNGB associated with SSB at the scale /, of the order of the seesaw 
scale. In this case the pNGB mass generated at one-loop can be written schematically as 
^™e ~ (y'^?/^t'^log(A^//i^)/(16'7r^). The mass mg is not proportional to /, but rather to the 
electroweak scale v. Besides, it can be much smaller than v because of the four powers of 
Yukawa couplings and the loop suppression. This is true even with a huge cutoff A, since the 
dependence on it is only logarithmic. However, at higher orders the quadratic divergence 
may reappear, e.g. through a two-loop diagram with a virtual Higgs exchanged across the 
neutrino loop, leading to m^ ~ g^y'^A'^/{lQ7c'^)'^. This result is not surprising, since the 
mechanism we adopted explains why the pNGB mass is related to the EW scale, but it does 
not address the stability of the EW scale against radiative corrections. In other words, as 
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already remarked by Hill and Ross ^, the Higgs sensitivity to quadratic corrections also 
enters in mg at higher order. In order for this two-loop correction to be negligible with respect 
to the one-loop estimate, one needs the cutoff of the Higgs boson loops, Ah, to be smaller than 
~ Anv. (By enlarging the global symmetry, it may be possible to remove also the two-loop 
quadratic divergence, see e.g. [16], but in general higher orders reintroduce the problem.) 
If one wants to stabilize a theory which includes a scale much larger than v (e.g. the scale 
/), one must address the usual hierarchy problem, e.g. postulating supersymmetry broken 
at Asusy ~ TeV, or a strongly interacting sector that condenses at Ac ~ TeV generating 
dynamically the EW scale. In this paper we assume the stability of the EW scale is realized, 
and thus we can adopt our one-loop estimates for the pNGB mass and couplings, in order 
to study phenomenology. 

4.1 A seesaw model leading to the Higgs portal 

Let us consider in some detail a specific U{l)x symmetry of the neutrino sector, such that the 
pNGB 6 acquires radiatively a coupling to the SM Higgs. Take the sterile neutrino charges 
X(z/^) = —1 and ^(z/g) = 1, and a scalar field $ with charge X{^) = 2, whose vev breaks 
the symmetry spontaneously. The interaction lagrangian involving the sterile neutrinos and 
the pNGB e reads 

Here M12 is a mass term allowed by U{l)x, while Mn and M22 are generated after SSB. The 
symmetry is broken explicitly by either ruai or ma2, depending on the U{l)x charge assigned 
to the lepton doublet /qJj Note that 9 can be rotated away by rephasing the neutrino fields 
if rriai = or 171^2 = 0. In this case the symmetry is restored and 6' is a true NGB. The 
dependence on 6 can be removed also when two independent entries of M are vanishing, 
therefore the symmetry breaking effects must vanish also in this limit. As we discussed in 
the previous section, this need of four different couplings to break the symmetry is the key 
to cancel quadratic divergencesjj 

The interactions in Eq. ([33]) generate an effective potential for the pNGB 9. We assumed 
that Mij and rriai are real and positive (see discussion below) and we performed the one-loop 
computation of the effective potential, which can be written as 

V^ff = ^eOH^H + 0{e') . (34) 



^ The choice of X-charges is made to reahzed a model of type (iii) , in accordance with the classification of 
the previous section. This choice is not unique. An equivalent possibility is a U{l)x symmetry with charges 
X{v '{) = —1, X{v2) — 1, X{la) = 0, broken spontaneously by a scalar field of charge X^cf)) — 1. Then in 



Eq. (331 one should replace ra^i — >■ niaie^^/^ and 771^2 — >■ '71^26"'^^-'^, while the entries Mn and M22 are 
independent from 9 and represent the source of explicit symmetry breaking. It is easy to check that the two 
realizations are equivalent, since the lagrangian is the same up to a rephasing of the neutrino fields. 

^Note that the field 9 can be rotated away from Eq. (33 1 by phase redefinitions, as long as the Uf entries 



that break explicitly U{l)x are set to zero. This is slightly different from the general case discussed below 
Eq. (29), where ny + 1 entries must vanish, for the trivial reason that here the same phase field appears 



(with opposite sign) in two independent entries. The only physically relevant fact still holds: the product of 
four independent entries is needed to generate a mass for 9. 



We find that the quadratically divergent contributions cancel exphcitly, as expected, and the 
logarithmically divergent ones give 



1 Mi2(Mn + M22) E« ^al"^a2 , A' 
47r2 /2 ^;2 °S^2 



A^^ ^^^"^^^^^]/^^^"^^"';Y"^"^ log%, (35) 



with a renormalization scale /z ~ Mjj and up to tiny corrections of higher order in mai/Mki- 
This type of models is particularly predictive, because the pNGB mass is generated by the 



same loops that generate A, that is, m^ is obtained from Eq. (34) by replacing the Higgs 
with its vev, rrig = Xv"^. By taking a common value mj\j = gf/\/2 [m = yv) for each entry of 
the Majorana (Dirac) neutrino mass matrix, Eq. (35) reduces to Eq. (|3]), up to a factor 2nf 



accounting for the sum over flavour indexes. As we saw in section [3} the value of the coupling 
A has a crucial role for the (partial) thermalization of 6 and thus for the determination of 
its relic density. 

A comment on the CP symmetry is in order. In general, the mass terms Mij and rriai may 
be complex, that is, they may carry phases that cannot be removed by a redefinition of the 
fields. These phases correspond to an explicit breaking of CP. In this case 6 would have both 
scalar and pseudo-scalar couplings to the fermions, with characteristic phenomenological 
signatures p]. Moreover, the effective potential would contain terms odd in 9, such as 
fiOWH (with fi '^ v), that induce a vev for 6 and a small mixing with the Higgs. This 
may endanger the stability of 6, since it would couple linearly to SM particles. Still, these 
couplings can be very small and in addition the 6 decays may be kinematically forbidden 
if me is sufficiently small. In this paper we do not investigate this more complicated CP- 
violating possibility, and we rather assume that CP is a good symmetry of the seesaw sector. 



In this case Mij and rriai in Eq. (33) are all real and 9 preserves its pseudo-scalar nature, 
since Eq. (33 ) is invariant under a CP transformation with 9 — t- —9. Usually the DM stability 
requires an additional (discrete) symmetry that forbids its decay into SM particles. In the 
present scenario the CP symmetry itself guarantees that 9 couples linearly only to the heavy 
sterile neutrinos, and we will see in section |5] that this makes 9 sufficiently long-lived. 

5 Constraints on the pNGB lifetime 

We first derive the couplings of 9 to the SM fermions and gauge bosons, and compare its decay 
width into light SM particles with the lifetime of the universe. Subsequently, assuming that 
^-particles account for the whole DM density, we discuss the more stringent astrophysical and 
cosmological constraints which exist on the decay of 9 to neutrinos, electrons and photons. 

5.1 9 couplings to the SM particles 

The pNGB lifetime is determined by its couplings to light fermions. These come from the 
9-N interaction, through the u-N mixing induced by the Dirac neutrino masses. In this way 
9 can decay to light neutrinos (at tree level) and to charged fermions (at one-loop level). In 
turn, these couplings to SM fermions could induce, through triangle loop diagrams, couplings 
to SM gauge bosons. These decays should be sufficiently slow to make the DM lifetime longer 
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than the age of the universe, tq ^ 5- lO^^s. This provides interesting constraints on the seesaw 
parameters, as we now discuss. 

Neutrinos. At tree level, 9 couples only to light neutrinos as follows: 

Ce.u = \Y. ^^^^75^^/3^ , Map = - 5^X,,(mM-i)„,M,,(M-im^),-;3 • (36) 

a/3 •' ij 

Here Xij is the power of e*^'-^ associated with the sterile neutrino mass matrix entry Mij. In 
the singlet Majoron model Xij = 1 for all i,j, therefore one finds fi,y = nii, = —mM'^rn^ , 
which is the usual seesaw formula. On the contrary, in our models based on a family- 
dependent U{l)x symmetry, 6 couples differently to each entry of M , as for example in the 



model of Eq. (33). This lagrangian leads to a total decay width (into both neutrinos and 
antineutrinos) 

T[e -^ vv) = Y^9euu^e , Qe^u = « , (37) 

where we neglected {m^i/nigy corrections in the phase space, mj,j being the light neu- 
trino mass eigenvalues. In the Majoron case one obtains gj^^j^ = Ylii^^til P- This width 
is smaller than I/tq for ge^y < 3 x 10~^^(MeV/m6i)^/^. Such tiny coupling is natural for 
6, because the pNGB couplings are suppressed by the SSB scale /: one needs / ^ 3 x 
10^GeV(m,y/eV)(rri5i/MeV)^/^, where we used the one-family approximation, gg^y ^ rn^lf. 
To translate this bound in the {me — g) plane we insert the relation / = ^J2vaj^l g in Eq. M|, 
which gives 
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(^) (l^)'(^)^- (-, 



geyy ~ 10" 

The condition l/r(^ — !■ vv') > tq excludes the region above the blue dashed line in Figs. 314 
Charged fermions. The pNGB 6 couples also to charged fermions, through EW one- 
loop diagrams. This effect arises because of the mixing between the sterile and the weakly- 
interacting neutrinos, in particular it is also operative in the singlet Majoron model [I] IT7]. 
The coupling to quarks is induced by a one-loop 6 — Z mixing diagram, with neutrinos in the 
loop. The coupling to charged leptons is generated by an analog Z-exchange diagram plus 
a triangle diagram with Pi^-exchange. The resulting coupling can be written in a compact 
form in the one-family approximation, as follows: 

r ■ -f fa ^ ^P V2mN , . 

C-eff = Wefffl5f e , geff = ±j—^^^mfm, , (39) 

where the sign is + (— ) for up quarks and charged leptons (down quarks), and Gp is the 
Fermi coupling constant. This effective coupling can be suppressed only taking a small 
g = \/2mN/ f, since all the other parameters in ggj-j can be determined experimentally. The 
decay width is given by 



i2 \ ^/2 
r(^ ^ //) = ^gtffme ( 1 - -/ . (40) 



1 2 _ A 4my 



57r •'•' \ m: 
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We remark that all models where the singlet Majoron is given a mass larger than one MeV 
and plays the role of DM candidate are constrained (or already excluded) by such decays 
into charged fermions. 

In our scenario, the requirement l/r(6' — > //) > tq leads to an upper bound g < 
2{MeV / mgY^'^ {eV / m,^) (MeV/mf), where we have assumed mg ^ 2mf for simplicity. This 
bound excludes the region above the red dashed line in Fig. |4j As a consequence, taking 
into account that g and mg are related by Eq. ([s]), and that the Dirac mass m = yv cannot 
be larger than ~ TeV, we conclude that mg close to the EW scale would imply a lifetime 
shorter than tq. In fact, even stronger constraints on the decay width into charged fermions 



come from astrophysical and cosmological observations (see section 5.2), which will lead to a 
stronger constraint, mg < 1 GeV. In particular, in this class of models 6 cannot be the ~ 50 
GeV DM candidate produced by the freeze-out of the A interaction with the Higgs, discussed 



in section 3^ One is left with the possibility of a sub-GeV DM candidate, because in this 
case the decays into charged fermion pairs are sufficiently slow (or forbidden kinematically) , 
and the correct relic density can be generated by the freeze-in mechanism. 

In the realistic three-family case, the coupling ggjj is generically of the same order, but 
with a complicated dependence on flavour parameters. In particular one may argue that 
some cancellation can take place, to reduce the 6 decay width. In addition, we have seen 
in section |4] that the pNGB mass generation depends crucially on the interplay between 
the flavour structures of the matrices m and M. Therefore, if one were to study the whole 
parameter space, one should know the explicit dependence of ggfj on the mass matrix entries. 



For illustration, we display the result for the model in Eq. (33), considering for simplicity 
only one lepton doublet. Writing the Dirac mass matrix as m = (r/ii 7712) and keeping terms 
up to order {mi/Mjk)"^, the effective coupling of 6 to fermions through the mixing with the 
Z gauge boson is given by 






^off = * ^^1 - -^F2 -i^^ nif T,fh,f e , (41) 



F, = m\ {ch^K + ^ + ^) - -^ c^s^K - m.m. cs ((c^ - s^)K - ^ + ^) , (42) 

F, = ml (cls^K + ^ + ^) - -? c^s^K + m,m, cs ((c^ - s^)K + ^ - ^) , (43) 

where T^f is the third isospin component of the left-handed part of the fermion /, Mi^2 are 
the eigenvalues of the matrix M, which is diagonalized by a rotation of angle S, deflned by 
tan 25 = 2Mi2/(Mii — M22), and we denoted c = cos 6 and s = sin 5. Finally, the loop 
function K = K{Mi, M2) is given by 

Mf + 4M1M2 + M^^ A{Ml + M1M2 + Ml) M, , . 

^ '' '^ ~ MiM2(Mi + M2) (Ml - M2)(Mi + M2Y ^ M2 ■ ^ ' 

Concerning the 9 decay into charged leptons, one needs to add the contribution of the W- 
exchange diagram. The main phenomenological constraint comes from 9 — >■ e+e~. Therefore, 
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in our simplified calculation we identify the active neutrino with the electron neutrino z/g. 
Then, the additional contribution to the effective coupling of 6 to electrons is given by 

^eee = ^ {-Y " ~T ) (47r)2 "^^^^se e , (45) 

that carries a relative factor —2 with respect to the Z-exchange contribution. This is why 



the sum of the two contributions in Eqs. (41) and (45) gives goee/fT^e = 98uu/fT^u = ~9eddl^^d 



consistently with the one-family result in Eq. (39). Also, the functions Fi 2 are of the order 
mj/Mj ~ nil,, so the couplings are of the same order as in Eq. (39). Still, cancellations 
between the various terms are possible leading to a suppression of the 6 decay width for 
special flavour structures. This possibility may deserve a future investigation, since an 
appropriate family symmetry could in principle raise the pNGB lifetime, so that 6 could 
become a viable DM candidate even for mg > 1 GeV, a region where freeze-out could 
lead to the observed relic density. In this paper we do not invoke a family symmetry for 
the suppression of Qq^j to happen, but rather we rely on the one-family estimate given in 



Eq. (jSgj). 

Gauge bosons. The last effective coupling of the pNGB that may have important 
phenomenological consequences is the one to photons. Let us discuss first the limit of exact 
global symmetry, and then comment on the effect of explicit breaking. In general, the 
couplings of a NGB to gauge bosons are controlled by the gauge anomalies of the associated 
global symmetry. In the case 6 is the Majoron, the global symmetry can be identified with 
B — L, because 6 can be rotated away from the Yukawa interactions by rephasing all the SM 
fermions with a. B — L transformation. Since B — L is anomaly-free with respect to the SM 
gauge symmetries, this rephasing does not generate anomalous couphngs of the type 9FF, 
rather the only leftover interaction is derivative, {0/ f)df^J'^_j;^, which is generated by the 
redefinition of the fermion kinetic terms. One concludes that the Majoron has no anomalous 
couplings to gauge bosonsjj 

This discussion is easily generalized to the U{l)x symmetry that we consider in this 
paper. First, note that the sterile neutrino charges are irrelevant for gauge anomalies. In 
order to allow for fermion masses, we can take X{la) = — X(e^) = Xa for a = e,fi,T, and 
zero X-charge for the quarks. Then it is easy to check that U{l)x is anomaly-free with 
respect to electromagnetism and colour, and thus 9 does not couple to photons and gluons. 
The anomaly with EW interactions is proportional to Yla ^" and can also be taken equal to 
zero for simplicity. In summary, in the limit of exact symmetry our NGB does not couple to 
gauge bosons. 

It is far more difficult to compute the pNGB couplings to gauge bosons in the presence 
of explicit symmetry breaking sources. One expects that such couplings will arise at some 
level, since there is no symmetry arguments to prevent them at all orders. However, the 



^ Equivalently, one may remove 9 from the Yukawa interactions by an L transformation; in this case one 
is left with the couphngs {9/ f)[d^J'l + {cwFwFw + cy Fy Fy) / {IQir"^)], where cw,y t^ account for the EW 
anomahes of the lepton number symmetry. There is no physical difference between the two pictures. In 
particular one can check that (i) 9 does not couple to weak gauge bosons at one-loop even in the L case, 
because the triangle diagrams generated by the derivative coupling to Jj^ cancel Cw,y exactly; (ii) 9 does 
not decay to two on-shell quarks even in the B — L case, because the quarks are vector-like under B ~ L. 
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Figure 3: The constraints on the DM lifetime in the {mg — g) plane, for nii, = 0.05 eV. The thick 
black curve, corresponding to the correct DM relic density, as well as the dotted lines, corresponding 
to constant values of rriN, were already presented in Fig. [!} The curve is dashed for mg < 1 keV, 
because DM is warm in this region (see the text), and for g < 10~^, because of the theoretical 
bound / < Mp, shown in Fig. m The blue solid (dashed) line is the upper bound on g from DM 
decays into neutrinos coming from astrophysics and cosmology (from the universe lifetime): the 
blue shaded region is therefore excluded. The red solid line is the analog bound for DM decays 
into e^e~: the red shaded region is correspondingly excluded. Finally, the brown dot-dashed line 
is a conservative estimate of the upper bound on g from DM decays into photons (see the text). 



computation of the lowest order non-vanishing contribution is non-trivial. We have shown 
explicitly above that, in the present framework, 6 couples at tree- level to neutrinos only, 
while, at one-loop, couplings to the Higgs boson as well as to charged fermions are induced. 
At two-loop there is a number of diagrams that connect 6 to two gauge bosons. These 
diagrams do not necessarily add up to zero: on the one hand, they involve the neutrino mass 
parameters that break explicitly U{l)x, thus they are not expected to respect the anomaly 
argument above; on the other hand, however, these parameters may not be sufficient to 
induce an operator 9FF already at two-loop order. 

The computation of these two-loop diagrams is beyond the purpose of the present paper, 
because we will show that a decay — )■ 77 induced at this order would be irrelevant for 
phenomenology anyway: once the constraints from ^ ff are imposed, the surviving DM 
parameter space is not further reduced by constraints on photons. To see this, we estimate 



the size of a two-loop contribution by taking the effective one- loop couplings in Eq. (39 ), and 



computing the usual fermion triangle diagrams with two final state photons. The effective 
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Figure 4: The same as in Fig. ^ but for a heavier neutrino mass scale, vtiy = 1 eV. The extra red 
dashed Hne is the upper bound from 1/V{6 — )• e^e~) > tq (in the case of Fig. ^the analog line lies 
entirely in the region 5 > 1 and mg > 1 GeV, which is not displayed). 



lagrangian can be written as 



£, 



077 



-^Qe-yyO^fiuafsF^'^F"'^ 



(46) 



and the corresponding decay width is 



T{e ^ 77) 



647r- 



^77 



TTlfi 
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y6»77 — 
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2m 



(47) 



arcsin x\/x [IH]. Note that even in the hmit of exact NGB, with me ^ m,f 
~ 1, the sum does not give zero when the sign in Eq. (39 ) is taken into account. 



where G{x) = 

and thus G{x) 

Therefore, we stress again that this is just a conservative estimate, and the computation 

of all the two-loop contributions may lead to a further cancellation. Neglecting the mass 

dependence in G, the requirement r(^ — >■ 77) < I/tq gives the order of magnitude constraint 

g < 4xlO^(MeV/?Tie)^/^(eV/mi,). This bound is weaker than the one obtained below Eq. (40) 



from r(^ — )■ //), as long as me < 1 GeV, which is the region relevant in the present scenario. 



5.2 Cosmological and astrophysical bounds on 9 couplings 

In practice, the DM lifetime has to be larger than the lifetime of the universe, because 
late DM decays affect several cosmological and astrophysical observations. To derive the 
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corresponding bounds on our scenario we make the assumption that 9 is all the DM in the 
universe. For the sake of simplicity, here too we take the approximation of one lepton family, 
barring large cancellations between the flavour parameters. A sub-GeV particle could decay 
into neutrinos, electrons and photons. We discuss these decay channels in turn. 

Let us first consider the decay 6 — )• pp. As it was noted in Ref. [5] such a decay could 
affect the expansion history of the universe, because it represents energy transfer from a non- 
relativistic to relativistic neutrinos. Using SNIa and CMB data, one obtains the bound 

m 

T{e -> z/z/) < 4.5 X 10^2° s-^ (or r > 700 Gyr) . (48) 



Using Eq. (37), the corresponding bound on the coupling of 9 to neutrinos is 



«..<4X10-(M5X)"\ (49) 

This upper bound can be improved in the range 30 MeV < mo < 200 MeV [20], using 
searches for the diffuse neutrino supernova background by Super-Kamiokande. For masses 
mg > 200 MeV the best limit come from atmospheric neutrino observations. Since the 
observed spectrum coincides with theoretical estimates, one can set an upper bound on gg^^ 
[201 121] • We use the relation in Eq. ( [38] ) to translate these constraints in the {mg — g) plane. 
The blue region shown in Figs. [Hpf is excluded by the combination of the observational 
bounds discussed above. 

Let us now consider the 6 — )■ e+e^ decay. We adapt to our case the analysis performed in 
Ref. [21]. For m,g ^ 20 MeV, the dominant physical process which constrains the parameter 
space is annihilation at rest, contributing to the 511 keV line. The limit is approximately 
given by 

(^^\ r(^ ^ g-g+) < 5 X 10-2^ s-^ . (50) 

\ mg J 



Using Eq. (40), this limit leads to 

goee < 9 X 10-2^ 1 ^ . (51) 

For 20 MeV ^ m,g ^ 1 GeV the dominant process is internal bremsstrahlung, i.e., photons 
radiated from the final electron or positron. One gets a slighly more stringent bound than 



follows from Eq. (39), is 



the one shown in Eq. (50) [21j. The predicted value of the coupling gg^^ in our model, as 



^P ™ ™ _ Q .. 1 n-23 ^ Q \ m,^ 



J(4j"'""'. = 3x10-"(^J^. (52) 



The corresponding observational bound is plotted in Figs. 3]|4 as the red shaded region. 



Note that the value m^ ~ 3 MeV, predicted by the Higgs portal, could lead to an excess 
of 511 keV 7 rays from the galactic center of the Milky Way. As discussed in Ref. [22], a 
~ 1.5 MeV electron or positron, produced by a 6* — )■ e+e" decay, is enough non-relativist ic 
to subsequently annihilate at rest in the galactic center, leading to a 7 line at 511 keV. In 
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particular, with such a mass one can easily obtain a 7 flux excess of the order of the one 
observed by the INTEGRAL 7 observatory [221 EH] . However, unless the DM galactic profile 
is much more cuspy than the usually considered profiles, a DM decay gives a flux that is not 
sufficiently peaked around the galactic center |22l Elj to be able to reproduce the morphology 
of the signal observed by INTEGRAL. 

We discuss now the decay 6 — )■ 77. For small masses, me < keV, the photons energy 
is absorbed by the baryonic gas in the early universe and the processes of recombination 
and reionization are affected. One can use the analysis done in Ref. [25] (see also pBJ) to 
bound the coupling ge-^-y- For larger masses, me > keV, photons are no longer absorbed and 
propagate freely. Their contribution to the isotropic diffuse photon background allows also 
to bound ge'y'y [23 [2Z] • 

For the region we are most interested in, stronger bounds can be obtained from the 
gamma-ray line emission limits from the Milky- Way central region. Indeed, in the range 40 
keV< me <16 MeV P3> one obtains 

r(«^„)<10--s-'(^). (53) 



Using Eq. (47), this corresponds to 



ge,, < 4 X 10-21 (^/^\ q^^-i (54) 

\ nig J 



Taking the estimate for ge-yy given in Eq. (47), all these constraints translate in an upper 



bound on g, which is shown in Figs. 3]|4 by the brown dot-dashed curve. 

We should mention that there are astrophysical constraints on the couplings geuv [21], 
geee [29J, and ge-j-y |2S1 ED], based on stellar energy loss due to 6 emission, provided m^ is 
low enough to be produced in stellar interiors. These limits are valid without the need to 
assume that 9 is DM. However, these astrophysical limits are in general weaker than the 
ones mentioned above. 

We also would like to note that the region we consider for pNGB DM masses includes 
the keV scale, corresponding to warm dark matter (WDM). Since at the epoch of structure 
formation WDM has free-streaming lengths below the Mpc scale, having WDM at least 
as a non-negligible DM component can alleviate some of the disagreement between the 
standard cold DM scenario and a variety of galactic observations at small scales ^31j|. The 
ideal observation which could place a limit on WDM is Lyman-a forest, i.e. the Lyman- 
a absorption produced by intervening neutral hydrogen in the spectra of light emitted by 
distant quasars. Using Ly-a observations together with other cosmological data sets, different 
groups have put lower limits on my/ dm-, assuming WDM has a thermal distribution and that 
it is the whole of DM. The bounds in the literature [32] differ by factors of a few, ranging 
from 0.5 keV up to 4 keV. In our case, the DM species is not in thermal equilibrium and the 
bounds can slightly change compared to the WDM thermal relic case. In addition, if WDM 
is only a part of the whole of DM the bounds are relaxed. We should mention that there 
are some potential problems in the obtention of these bounds, and they should be regarded 
as controversial. On the one hand there could be large systematic errors, and on the other 
hand the Ly-a analysis has to be performed at scales which are already in the non-linear 
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regime, where calculations are less reliable. If one disposes of the Ly-a data these bounds 
on WDM disappear altogether [33] . 

6 Conclusions 

We proposed a new pseudo-scalar gauge-singlet DM candidate 9 with mass in the keV - 
MeV range. Its couplings to the SM particles are feeble, because they are mediated by new 
physics at a large scale /, which we identify with the seesaw scale. The 9 relic density can be 
produced by the scattering with the heavy particles (the sterile neutrinos), at temperatures 
of the order of /, or alternatively it can be generated at the EW scale through the Higgs 
portal, by the tiny 9-H coupling A ~ 10~^°, which is induced by the seesaw interactions. 
Today, 9 decays into light neutrinos and, if it is heavier than one MeV, into e^e~ , with rates 
that can saturate the present upper bounds. 

We argued that such a candidate is theoretically well-motivated. The heavy new physics 
sector is generically associated with several global symmetries. Some of the corresponding 
pNGBs may remain light, if the explicit symmetry breaking effects are sufficiently small. For 
concreteness, we demonstrated that a U{l)x family symmetry of the sterile neutrino sector 
can be broken collectively by a set of neutrino Yukawa couplings, so that the pNGB mass 
mg is proportional to the EW scale and to the (small) product of the Yukawa couplings. In 
such scenario tuq is not quadratically sensitive to the cutoff (at leading order), therefore the 
presence of a light DM scalar below the EW scale is justified. 

In order to calculate the 9 relic density, we computed the rate of the 9 interactions with 
a sterile neutrino N and, through the Higgs portal, with the SM particles, and we studied 
numerically the Boltzmann equation for the 9 number density. In this framework there 
are only two independent parameters, the mass of the pNGB, m^, and the coupling of the 
pNGB to the sterile neutrino, g, a feature which makes our scenario especially predictive. 
We find that the Higgs portal produces the desired DM relic density, through the freeze-in 
mechanism, for a unique value of the DM mass, mg ~ 3 MeV. This prediction relies only 
on the relation between the DM mass and its coupling to the Higgs, m^ = Af^, which is 
well-justified in the case of our pNGB. As long as the reheating temperature is sufficiently 
high, the sterile neutrino portal can produce the required relic density by freeze-in for smaller 
values of mg, from ~ 3 MeV down to 0.15 keV. These results are summarized in Fig. [1} 

The constraints from r(^ — ?■ vv), however, exclude the region mg < 1(100) keV, for a 
light neutrino mass rrti, = 0.05(1) eV. In addition, the constraints from r{9 — )■ e~^e~) put 
an upper bound on the coupling g in the region mg > 1 MeV: for rriu = 0.05(1) eV and 
g ~ 10^'^ (lO"'*), the expected electron-positron fiux is close to the present sensitivity. These 
results are summarized in Figs. 3]|4 



In turn, the seesaw scale can be constrained, by requiring our pNGB to be a viable DM 
candidate in agreement with all the bounds above. In the case of freeze-in through the sterile 
neutrino portal, we find 10^ GeV < m^v ^ 10^° GeV, and a SSB scale / ~ lO'^m^v. In the 
case of the Higgs portal, one has instead 10^° GeV < m^v < 10^^ GeV, and correspondingly 
10^3 GeV <f<Mp. 

Finally, we note that our analysis of the parameter space was mostly performed in the 
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one-family approximation, that is, neglecting possible hierarchies among the sterile neutrino 
mass parameters and the neutrino Dirac Yukawa couplings. Some of our main results are 
independent from this approximation, such as the allowed range for the DM mass. On the 
contrary, the bounds on g and itim are obviously sensitive to the flavour structure. A rough 
idea of these dependence can be grasped by comparing the case rriy = 0.05 eV (Fig. [s]) with 
the case m,^ = 1 eV (Fig. 111). A more detailed exploration of the flavour parameter space of 
this scenario will be desirable, in particular, if one wants to compare with neutrino flavour 
models. 

Acknowledgements 

It is a pleasure to thank K.S.Babu, X.Chu, J.R.Espinosa, D.Greynat, K.Jedamzik, S.Lavi- 
gnac, S. Peris, S.Palomares-Ruiz, A.Pomarol, M.Quiros, J.Redondo, J.Serra, M.H.G.Tytgat 
and A.Varagnolo for very useful discussions. The work of MF was supported in part by the 
Spanish Research Project CICYT-FEDER-FPA-2008-01430, the Consolider-Ingenio 2010 
Programme CPAN (CS-D2007-00042), and the Marie-Curie Reintegration Grant PERG06- 
GA-2009-256374 within the European Community FP7. MF also thanks the Service de 
Physique Theorique of the Universite Libre de Bruxelles for hospitality. The work of TH is 
supported by the FNRS-FRS, the IISN and the Belgian Science Pohcy (lAP VI-11). TH 
thanks J.R.Espinosa and the IFAE at Barcelona for hospitality. The work of EM has been 
partly supported by CICYT-FEDER-FPA2008-01430 and by 2009SGR894. 

References 

[1] Y. Chikashige, R. N. Mohapatra and R. D. Peccei, Phys. Lett. B 98 (1981) 265. 

[2] E. K. Akhmedov, Z. G. Berezhiani, G. Senjanovic and Z. j. Tao, Phys. Rev. D 47 (1993) 
3245 tarXiv:hep-ph/9208230i . 

[3] L Z. Rothst ein, K. S. Babu and D. Seckel, Nucl. Phys. B 403 (1993) 725 |arXiv:hep- 
|ph/93012l3l . 



[4] V. Berezinsky and J. W. F. Valle, Phys. Lett. B 318 (1993) 360 |arXiv:hep-ph/9309214| . 

[5] M. Lattanzi and J. W. F. Valle, Phys. Rev. Lett. 99, 121301 (2007) |arXiv:0705.2406l 

[astro-ph]]. 

[6] F. Bazzocchi, M. Lattanzi, S. Riemer-Sorensen and J. W. F. Valle, JCAP 0808 (2008) 
013 [arXiv:080 5.2372 [astro-ph]]; J. N. Esteves, F. R. Joaquim, A. S. Jos hipura, J. C. Ro- 



mao, M. A. Tortola and J. W. F. Valle, Phys. Rev. D 82 (2010) 073008 |arXiv: 1007.0898 
[hep-ph]]. 

[7] P. H. Gu, E. Ma and U. Sarkar, Phys. Lett. B 690 (2010) 145 f arXiv:1004.19l9) [hep-ph]]. 

[8] L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 1003, 080 (2010) 
[arXiv:0911.1120 [hep-ph]]. 

28 



[9] C. T. Hill and G. G. Ross, Phys. Lett. B 203 (1988) 125; C. T. Hill and G. G. Ross, 
Nucl. Phys. B 311 (1988) 253. 



[10] M. Schmaltz and D. Tucker-Smith, Ann. Rev. Nucl. Part. Sci. 55 (2005) 229 |'arXiv:hep- 
ph/0502182|; M. Perelstein, Prog. Part. Nucl. Phys. 58 (2007) 247 |ar Xiv:hep- 

[ph/0512128|. 



[11] W. L. Guo and Y. L. Wu, JHEP 1010 (2010) 083 [arXiv: 1006.2518 [hep-ph]]. 



[12] M. Farina, D. Pappadopulo and A. Strumia, Phys. Lett. B 688, 329 (2010) 



tarXiv:0912:5038| [hep-ph]] 



[13] C. E. Yaguna, JHEP 1108 (2011) 060 |arXiv:1105.1654 [hep-ph]]. 



[14] A. Y. .Smirnov and R. Zukanovich Funchal, Phys. Rev. D74, 013001 (2006) [ hep- 
fph/0603009]. 



[15] P. Sikivie, Lect. Notes Phys. 741 (2008) 19 [arXiv: astro-ph/ 061 0440 



[16] R. Bar bieri, L. J. Hall, S. J. Oliver and A. Strumia, Phys. Lett. B 625 (2005) 189 
|arXiv:hep-ph/0505124l . 

[17] A. Pilaftsis, Phys. Rev. D 49 (1994) 2398 |arXiv:hep-ph 79308258] . 



[18] L. AmetUer, L. Bergstrom, A. Bramon and E. Masso, Nucl. Phys. B 228 (1983) 301. 
[19] Y. Gong and X. Chen, Phys. Rev. D 77, 103511 (2008) [arXiv:0 802.2"296 [astro-ph]]. 
[20] S. Palomares-Ruiz, Phys. Lett. B 665, 50 (2008) [arXiv: 0712. 19 37 [astro-ph]]. 



[21] N. F. Bell, A. J. Galea and K. Petraki, Phys. Rev. D 82, 023514 (2010) [ arXiv: 1004. 1008 

[astro-ph.HE]]. 



[22] N. Prantzos et ai, arXiv: 1009.4620 [astro-ph.HE], and Refs. therein 



[23] G. Weidenspointner, G. Skinner, P. Jean, J. Knodlseder, P. von Ballmoos, G. Bignami, 
R. Diehl, A. W. Strong et al, Nature 451 (2008) 159-162. 

[24] Y. Ascasibar, P. Jean, C. Boehm and J. Knoedlseder, Mon. Not. Roy. Astron. Soc. 368 

(2006) 1695 |arXiv:astro-ph/0507142| . 

[25] L. Zhang, X. Chen, M. Kamionkowski, Z. g. Si and Z. Zheng, Phys. Rev. D 76, 061301 

(2007) |arXiv:0704.2444 [astro-ph]]. 



[26] X. L. Chen and M. Kamionkowski, Phys. Rev. D 70, 043502 (2004) | arXiv:astro- 
lph/0310473]. 



[27] H. Yuksel and M. D. Kistler, Phys. Rev. D 78, 023502 (2008) |arX iv:0711.2906| [astro- 

ph]]. 



29 



[28] J. A. Grifols, E. Masso and S. Peris, Phys. Lett. B 215 (1988) 593. 

[29] G. G. Raffelt, "Stars As Laboratories For Fundamental Physics: The Astrophysics Of 
Neutrinos, Axions, And Other Weakly Interacting Particles," Chicago, USA: Univ. Pr. 
(1996) 664 P 



[30] E. Masso and R. Toldra, Phys. Rev. D 52 (1995) 1755 | arXiv:hep-ph/9503293 

[31] B. Moore, S. Ghigna, F. Governato, G. Lake, T. R. Quinn, J. Stadel and P. Tozzi, As- 
trophys. J. 524, L19 (1999); G. Genti le, A. Burkert, P . Salucci, U. Klein and F. Walter, 
Astrophys. J. 634, L145 (2005) tarXiv:astro-ph/0510 60r| . 

[32] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese and A. Riotto, Phys. Rev. D 
71, 063534 (2005) |arX iv:astro-ph/0501562|; K. Abazajian, Phys. Rev. D 73, 063513 
(2006) |arXiv:astro-ph/0512631|; M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese 
and A. Riotto, Phys. Rev. Lett. 97, 071301 (2006) |arXiv:astro-ph/0605706|; U. Seljak, 
A. Makarov, P. McDonald and H. Trac, Phys. Rev. Lett. 97, 191303 (2006) [arXiv: astro-] 



ph/0602430| ; M. Viel, G. D. Becker, J. S. Bolton, M. G. Ha ehnelt, M . Ranch and 



W. L. W. Sargent, Phys. Rev. Lett. 100, 041304 (2008) |arXiv:0709.0131 [astro-ph]]. 



[33] S. Hannestad and G. Raffelt, JCAP 0404, 008 (2004) |arXiv:hep-ph/0312154| . 



30 



